;

; Open file.
file_name='G:\CALIPSO\PRO_CODE\data\CAL_LID_L2_VFM-Standard-V4-20.2019-07-12T04-47-04ZD.hdf'

fid = HDF_SD_START(file_name, /READ)

; Define data name. Use "HDFView" or "hdp" tool to check datasets
; available inside the file.

data_name="Feature_Classification_Flags"
index=HDF_SD_NAMETOINDEX(fid, data_name)

; Retrieve FCF data.
sds_id=HDF_SD_SELECT(fid, index)
HDF_SD_GETDATA, sds_id, fcf

; Retrieve lat data.
data_name="Latitude"
index=HDF_SD_NAMETOINDEX(fid, data_name)
sds_id=HDF_SD_SELECT(fid, index)
HDF_SD_GETDATA, sds_id, lat

; Retrieve lon data.
data_name="Longitude"
index=HDF_SD_NAMETOINDEX(fid, data_name)
sds_id=HDF_SD_SELECT(fid, index)
HDF_SD_GETDATA, sds_id, lon

; End access to SDS.
HDF_SD_ENDACCESS, sds_id

; Close file.
HDF_SD_END, fid

; Users need to understand the layout of the Feature_Classification_Flag
; dataset.
;
; The Feature_Classification_Flag values are stored as a sequence of 5515
;  element arrays (i.e., as an N x 5515 matrix, where N is the number of
; separate records in the file). In this file, N is 4224.
;
;  Each array represents a 5 km "chunk" of data,
; and each array element contains the feature classification information for a
;  single range resolution element in the Level 0 lidar data downlinked from
; the satellite. As shown in the summary below, the vertical and horizontal
; resolution of the CALIPSO data varies as a function of altitude above mean
; sea level (see Hunt et al., 2009).
;
; Here's the summary of number of prfofiles per 5 km.
;
; 3 profiles for 20.2km (base) to 30.1km (top) @ 180m
; (index 1-165 / 55 samples per profile)
; 5 profiles for 8.2km (base) to 20.2km (top) @ 60m
; (index 166 - 1165 / 200 samples per profile)
; 15 profiles for -0.5km (base) to 8.2km (top) @ 30m
; (index 1166 - 5515 / 290 samples per profile)
;
; 3 profiles mean horizontal resolution is 1667m because 3 * 1667m = 5km.
; 55 samples mean vertical resolution is 180m because 55 * 180m = 9.9km  =
; 30.1km - 20.2km.
;
; 1113.2m equals to 0.01 degree difference.
; 111.32m equals to 0.001 degree difference.
;
; Thus, we can ignore horizontal resolution for this global plot example.
;
; In summary, profile size determines horizontal resolution and sample size
; determines the vertical resolution.
;
; Each vertical feature mask record is a 16 bit integer.  See [1] for details.
; Bits | Description
; ----------------
; 1-3  | Feature Type
; 4-5  | Feature Type QA
; ...   |...
; 14-16 | Horizontal averaging
;
; In this example, we'll focus only on "Featrue type."
;
; However, the resolution of the height will be different.
;
; Altitude Lidar data is in "metadta" [2] stored as HDF4 Vdata.
; Lidar_Data_Altitudes has 583 records it does not match dataset size
; 565(=55+200+290).
; There are 5 below -0.5km and 30 above 30.1km.
;
; Therefore, we cannot not rely on the Vdata for altitude.
; Instead, we should calculate altitude from the data specification.
;
; For each 5515 at a specific lat/lon, we can construct cloud bit vector over
; altitude.
;
; For example, Feature_Classification_Flags[loc][55] means,
; Longitude[loc] and altitude = 30.1km.
;
; For another example, Feature_Classification_Flags[loc][56] means,
; Longitude[loc] + 1667m and altitude = 20.2km.
;
; There are many possibilites to plot this data.
; Here, we'll pick a specific altitude and plot Feature Type on
; 2-D map.
;
; Subset data at 2500m (= -0.5km + 30m * 100) altitude.
alt_index = 1256
data = fcf(alt_index, *)

; Select the first 3 bit for the feature type.
data = data AND '0007'X

; lat/lon are 2-D, Nx1 array.
lat = reform(lat)
lon = reform(lon)
data = reform(data)
dim=size(data,/dim)
longname = 'Feature Type at altitude=2500m'

; Define the color table.
levels = 8
;ct = COLORTABLE([[150, 0,   255,  0,   255,  200, 100, 50], $
;                 [150, 0,   255,  255, 0,    100, 50,  25], $
;                 [150, 255, 0,    0,   0,    255, 255, 125]], $
;                 NCOLORS = levels, /TRANSPOSE)
ct = COLORTABLE([[255, 224, 135, 255, 255, 46,  152, 0], $
                 [255, 255, 206, 215,  0,  139, 251, 0], $
                 [255, 255, 235,  0,   0,  87,  152, 0]], $
                 NCOLORS = levels,   /TRANSPOSE)                 

; Generate the plot.
m = MAP('Geographic', TITLE=file_name, FONT_SIZE=9, /BUFFER)
t1 = TEXT(0.35, 0.2, longname)

; MAGNITUDE requires BYTE type data.
data = BYTE(data)

; We use SCATTERPLOT because data is 2-d lat/lon swath.
c1 = SCATTERPLOT(lon, lat, OVERPLOT=m, $
                 MAGNITUDE=data, $
                 RGB_TABLE=ct, $
                 POSITION=[0.1, 0.1, 0.83, 0.9],$
                 /SYM_FILLED, SYMBOL='o', SYM_SIZE=0.1)
mc = MAPCONTINENTS()

; We need a custom colorbar because we use SCATTERPLOT().
cb = COLORBAR(RGB_TABLE=ct, RANGE=[0,8], ORIENTATION=1, BORDER=1,$
              TICKVALUES=[0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5], $
              TICKNAME=['invalid', 'clear', 'cloud', 'aerosol', $
              'strato.', 'surface', 'subsurf.', 'no signal'], $
              TEXTPOS=1, POSITION=[0.85,0.2,0.87,0.8], TITLE=units)

output_dir = 'G:\CALIPSO\PRO_CODE\data\'
base_name = file_basename(file_name)
file_core = strmid(base_name, 0, strlen(base_name) - 4)
output_file = output_dir + 'geo_wrold_track_' + file_core + '.png'
;png = file_name + '.idl.png'
c1.save, output_file, HEIGHT=600, WIDTH=800
;EXIT
end
; References
; [1] http://www-calipso.larc.nasa.gov/resources/calipso_users_guide/data_summaries/vfm/index.php
; [2] http://www-calipso.larc.nasa.gov/resources/calipso_users_guide/data_summaries/vfm/index.php#heading03

